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ABSTRACT 


The inner solar system's modern orbital architecture provides inferences into the epoch of terrestrial 
planet formation; a ~100 Myr time period of planet growth via collisions with planetesimals and 
other proto-planets. While classic numerical simulations of this scenario adequately reproduced the 
correct number of terrestrial worlds, their semi-major axes and approximate formation timescales, they 
struggled to replicate the Earth-Mars and Venus-Mercury mass ratios (79 and 15, respectively). In a 
series of past independent investigations, we demonstrated that Mars' mass is possibly the result of 
Jupiter and Saturn's early orbital evolution, while Mercury's diminutive size might be the consequence 
of a primordial mass deficit in the region (potentially the result of the growing Earth's early outward 
migration). Here, we combine these ideas in a single modeled scenario designed to simultaneously 
reproduce the formation of all four terrestrial planets and the modern orbits of the giant planets in 
broad strokes. By evaluating our Mercury analogs’ core mass fractions, masses, and orbital offsets from 
Venus, we favor a scenario where Mercury forms through a series of violent erosive collisions between 
a number of ~Mercury-mass embryos in the inner part of the terrestrial disk. We also compare 
cases where the gas giants begin the simulation locked in a compact 3:2 resonant configuration to 
a more relaxed 2:1 orientation and find the former to be more successful. In 2:1 cases, the entire 
Mercury-forming region is often depleted due to strong sweeping secular resonances that also tend 
to overly excite the orbits of Earth and Venus as they grow. While our model is quite successful 
at replicating Mercury's massive core and dynamically isolated orbit, the planets’ low mass remains 
extremely challenging to match. Indeed, the majority of our Mercury analogs have masses that are 2-4 
times that of the real planet. Finally, we discuss the merits and drawbacks of alternative evolutionary 
scenarios and initial disk conditions (specifically a narrow annulus of material between 0.7-1.0 au). We 
argue that the results of our N-body accretion models are not sufficient to break degeneracies between 
these different models, and implore future studies to apply further cosmochemical and dynamical 
constraints on terrestrial planet formation models. 


1. INTRODUCTION 


Models of the late stages of terrestrial planet for- 
mation (e.g.: Wetherill 1980a; Raymond et al. 2020) 
typically follow the dynamical evolution of a collec- 
tion of ~10-100 proto-planets (dubbed “embryos”, e.g.: 
Wetherill & Stewart 1993; Kokubo & Ida 1996, 2000; 
Morishima et al. 2010) engulfed within a distribution of 
~100-1,000 km planetesimals (e.g.: Youdin & Goodman 
2005; Johansen et al. 2015; Drazkowska & Dullemond 
2018; Lichtenberg et al. 2021). As the giant planets' 
gaseous compositions indicate that they formed within 
the lifetime of the Sun's primordial gas disk (1-5 Myr: 
Haisch et al. 2001; Hernández et al. 2007), their presence 
in simulations of the terrestrial world's ultimate accre- 
tion is essential (Chambers & Wetherill 1998; Levison & 


Agnor 2003). If the giant planets’ orbits remain circu- 
lar through the duration of terrestrial planet formation 
as predicted in hydrodynamical disk models (Papaloizou 
& Larwood 2000; Masset & Snellgrove 2001; Morbidelli 
& Crida 2007; Pierens & Nelson 2008; Zhang & Zhou 
2010), the final masses of Mars and the asteroid belt are 
too massive by at least an order of magnitude (Cham- 
bers 2001; Raymond et al. 2006, 2009; Lykawka & Ito 
2019; Woo et al. 2022). However, the moderate degree 
of radial mixing in such a scenario also has the advan- 
tage of aiding in the replication of Earth's water content 
(Raymond et al. 2004) and disparities between the iso- 
topic compositions of Earth and Mars (Tang & Dauphas 
2014; Dauphas 2017; Woo et al. 2021b). 

It is possible to generate a small Mars and low-mass 
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asteroid belt with circular giant planet orbits if, rather 
than spanning the full radial range between Mercury and 
Jupiter’s modern orbits, the majority of the terrestrial 
disk’s mass is concentrated in a narrow annulus between 
the current orbits of Venus and Earth (Agnor et al. 1999; 
Hansen 2009; Izidoro et al. 2014). Multiple explanations 
(see section 2.1) for these specific initial conditions have 
been proposed and robustly tested in the recent litera- 
ture (for a more complete summary see: Clement et al. 
2018; Raymond et al. 2020). Among others, these in- 
clude material removal during the gas disk phase via 
sweeping secular resonances with eccentric giant plan- 
ets (the “Dynamical Shake-up:” Nagasawa et al. 2000; 
Thommes et al. 2008; Bromley & Kenyon 2017; Woo 
et al. 2021a), Jupiter directly sculpting the terrestrial 
disk by migrating inward to ~1.5 au and back out due 
to interactions with the nebular gas (the “Grand Tack:” 
Walsh et al. 2011; Pierens & Raymond 2011; Jacobson 
& Morbidelli 2014; Brasser et al. 2016; Deienno et al. 
2016; Walsh & Levison 2016), or highly localized plan- 
etesimal formation (the “low-mass asteroid belt” or “de- 
pleted disk”: Izidoro et al. 2015, 2016; Drazkowska et al. 
2016; Raymond & Izidoro 2017b; Lykawka & Ito 2019; 
Mah & Brasser 2021; Morbidelli et al. 2022; Izidoro et al. 
2021b). 

An additional complication on this series of events is 
the giant planets’ acquisition of their modern, moder- 
ately eccentric orbits (Hahn & Malhotra 1999; Tsiga- 
nis et al. 2005) through an epoch of mutual encounters 
(Morbidelli et al. 2009; Nesvorny 2011). While a low- 
mass Mars is a regular outcome in simulations where the 
giant planets’ inhabit their current dynamical configura- 
tion for the duration of the simulation (Raymond et al. 
2009; Kaib & Cowan 2015; Lykawka & Ito 2019; Woo 
et al. 2021a; Nesvorny et al. 2021), such a scenario con- 
flicts with the predictions of many disk models of giant 
planet formation and early evolution (Morbidelli et al. 
2007; Pierens & Nelson 2008; Zhang & Zhou 2010), and 
also cannot explain Earth and Mars’ disparate compo- 
sitions (Woo et al. 2022). 

In the classic Nice Model (Gomes et al. 2005; Tsiganis 
et al. 2005; Levison et al. 2011), the giant planets’ orbits 
transition from circular to eccentric when they undergo 
an epoch of dynamical instability at t ~ 650 Myr; coinci- 
dent with a perceived spike in lunar cratering known as 
the Late Heavy Bombardment (Tera et al. 1974). How- 
ever, simulations of a late instability typically result in 
the catastrophic disruption of the fully formed terrestrial 
system (Brasser et al. 2009; Agnor & Lin 2012; Brasser 
et al. 2013; Kaib & Chambers 2016). Moreover, multi- 
ple recently derived geophysical (e.g.: Evans et al. 2018; 
Morbidelli et al. 2018; Mojzsis et al. 2019; Brasser et al. 
2020), geochemical (e.g.: Boehnke & Harrison 2016; Zell- 
ner 2017; Goodrich et al. 2021; Worsham & Kleine 2021) 


and dynamical (e.g.: Nesvorny et al. 2018; Quarles & 
Kaib 2019; Ribeiro et al. 2020; Nesvorny 2021; Liu et al. 
2022) constraints have been interpreted to strongly sug- 
gest that the instability happened within the first 100 
Myr after the solar system’s birth, and perhaps much 
earlier. 

While an instability at t ~ 10-100 Myr might have 
destabilized a nearly-formed terrestrial system (DeSouza 
et al. 2021) and triggered the giant impact that formed 
the Moon (Benz et al. 1986; Canup 2004b, 2012; Cuk 
& Stewart 2012) around the same epoch (730-100 Myr 
after nebular dispersal as inferred via isotopic dating: 
Wood & Halliday 2005; Kleine et al. 2009; Rudge et al. 
2010; Kleine & Walker 2017), an instability occurring 
before embryos at 1.5 au grow beyond a Mars-mass has 
been shown to reduce the final masses of Mars analogs 
and the asteroid belt (Clement et al. 2018; Deienno et al. 
2018; Clement et al. 2019c, 2021b; Nesvorny et al. 2021). 
In this paper, we continue to develop this model with 
new numerical simulations. In contrast to our past work 
(described below), our current effort specifically incor- 
porates reduced integration timesteps and inner terres- 
trial disk structures necessary for generating Mercury- 
like planets (Clement & Chambers 2021; Clement et al. 
2021a,c). While the distributions of planetesimals and 
embryos in the Mercury-region used here were found to 
be successful at producing reasonable Mercury analogs 
in a pervious series of studies (described below), it is 
currently unclear whether strong sweeping secular reso- 
nances with Jupiter in the inner solar system that occur 
during the giant planet instability would adversely affect 
Mercury's growth in such a scenario. For this reason, the 
primary goal of this paper is to understand how the Nice 
Model instability might have affected the formation of 
Mercury. 


2. BACKGROUND 


2.1. Models replicating Mars? mass 
2.1.1. Giant Planet Migration: The Grand Tack 


Perhaps the most recognizable terrestrial planet for- 
mation scenario, the “Grand Tack" model (as in the sail- 
ing maneuver: Walsh et al. 2011; Jacobson & Morbidelli 
2014; Brasser et al. 2016) supposes that the inward- 
outward migration (Pierens & Raymond 2011) of Jupiter 
and Saturn during the gas disk phase truncated the ter- 
restrial disk of planetesimals around ~1.0 au (Wether- 
ill 1978; Morishima et al. 2008; Hansen 2009). Specifi- 
cally, Jupiter's presence in the terrestrial region serves to 
both evacuate a large fraction of material from the Mars- 
forming region and simultaneously transfer volatile-rich 
C-type asteroids from the outer solar system to the as- 
teroid belt (DeMeo & Carry 2013). In spite of these 
consistencies, the mechanism for Jupiter’s tack strongly 
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Figure 1. Cartoon representation of the various initial conditions tested here, and in our past investigations of the early 
instability scenario in Paper 1, Paper 2 and Paper 3. The cartoon images of the outer planets are from www.kissclipart.com 


depends on unconstrained disk parameters, and has yet 
to be validated in simulations incorporating gas accre- 
tion (Raymond & Morbidelli 2014). Additionally, the 
strong radial mixing of material that occurs in the sce- 
nario is potentially inconsistent (Mah & Brasser 2021) 
with the disperate isotopic compositions of Earth (Lod- 
ders 2000; Dauphas 2017) and Mars (Tang & Dauphas 
2014). 


2.1.2. Secular Resonance Sweeping: The Dynamical 
Shake-up 


If Jupiter's orbit was eccentric while gas persisted 
in the solar nebula, the powerful v5 resonance would 
have swept from the belt's outer edge, inward to the 
vicinity Earth's orbit as the disk photo-evaporated (the 
^Dynamical Shake-up" described in: Nagasawa et al. 
2000; Thommes et al. 2008; Bromley & Kenyon 2017). 
While early hydrodynamical models (Masset & Snell- 
grove 2001; Morbidelli et al. 2007) indicated that such 
primordial eccentricities are unlikely outcomes of the gas 
disk phase, recent work suggests that primordial excita- 
tion is plausible for certain disk parameters (Kley et al. 
2004; Zhang & Zhou 2010; Pierens et al. 2014). More- 
over, the giant planets’ modern configuration (Clement 
et al. 2021d) and the disparate accretion zones of Earth 
and Mars (Woo et al. 2021b) are broadly consistent with 
such a scenario. 


2.1.3. Highly Localized Planetesimal Formation Efficiency: 
The Low-Mass Asteroid Belt 

Modern models of planetesimal formation have 
demonstrated how the process is highly sensitive to 
the local thermal and structural properties in the disk 
(Simon et al. 2016; Drazkowska et al. 2016; Lichten- 
berg et al. 2021). Moreover, ALMA observations show- 
ing non-uniform radial concentrations of dust in proto- 
planetary disks seem to support highly localized, or ra- 
dially dependent planetesimal formation. The low-mass 
asteroid belt model (Izidoro et al. 2015; Raymond & 
Izidoro 2017a,b; Izidoro et al. 2021a) proposes that very 
little solid material ever existed in the Mars-forming and 
primordial asteroid belt regions, and the terrestrial plan- 
ets grew from a correspondingly steep surface density 
profile of material. However, it remains difficult to de- 
termine how realistic such a distribution of solid mate- 
rial is. In this manner, it is challenging for disk models 
to demonstrate the contemporary formation of two iso- 
topically distinct populations of iron-meteorite parent 
body planetesimals at disparate radial locations (see: 
Lichtenberg et al. 2021; Morbidelli et al. 2022; Izidoro 
et al. 2021b; Chambers 2022). 


2.1.4. Terrestrial migration: Convergent and outward 


It is also possible that the terrestrial embryos them- 
selves migrated from their initial formation locations. 
Indeed, depending on the particular physical structure 
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and thermal profile of the nebular disk’s inner com- 
ponent, inward, 2Mars-mass embryos can experience 
inward, outward or convergent migration (e.g.: Cress- 
well & Nelson 2008; Lyra et al. 2010; Paardekooper 
et al. 2011; Bitsch et al. 2015; Eklund & Masset 2017). 
Recently, Broz et al. (2021) leveraged this concept to 
demonstrate that the masses and orbits of all four ter- 
restrial planets might be a consequence of a more diffuse 
collection of terrestrial embryos (~0.4-1.8 au) being re- 
shaped via convergent migration towards ~1 au. Simi- 
larly, Clement et al. (2021c) invoked outward migration 
of Earth and Venus’ precursor embryos to reconcile iso- 
topic differences between the Earth and Mars (Tang & 
Dauphas 2014; Dauphas 2017), as well as provide and 
explanation for Mercury’s diminutive mass and isolated 
orbit (Clement & Chambers 2021). However, these mod- 
els rely heavily on a priori assumptions of the solar neb- 
ula’s structure. 


2.2. The Early Instability Scenario 


In Clement et al. (2018, hereafter Paper 1) we studied 
the effects of the Nice Model instability on the forming 
terrestrial planets by varying the time at which the in- 
stability transpires within the accretion process. Here, 
and throughout this manuscript, we refer to the “insta- 
bility delay" as the amount of time a terrestrial planet 
formation simulation progresses before the instability 
ensues. We loosely correlate the beginning of the ac- 
cretion simulations with nebular dissipation (or 2-3 Myr 
after the formation of Calcium Aluminum-Rich Inclu- 
sions: CAIs); however, this connection is not exact. 

In all instability models tested, the terrestrial planets 
formed from an extended disk of 100 embryos and 1,000 
planetesimals with an outer boundary at 4.0 au (mod- 
eled after the classic initial conditions of Chambers & 
Wetherill 1998, hereafter C& W98, see figure 1), and the 
giant planet instability was modeled using the preferred 
initial conditions of Nesvorny & Morbidelli (2012, here- 
after referred to as the N&M12 model). Specifically, 
the N&M12 instability model assumes that Jupiter and 
Saturn emerged from the gas disk locked in a 3:2 mean 
motion resonance (MMR). In Paper 1, we considered 
any simulation where Jupiter and Saturn's final orbital 
period ratio (Ps/ P; = 2.49 in the modern solar system) 
finished less than 2.8 to be successful. We tested five 
and six planet giant planet models of this type (specif- 
ically, resonant chains of the form 3:2,3:2,3:2,3:2 and 


3:2,4:3,3:2,3:2,3:2'), and instability delays of 0.01, 0.1, 
1.0 and 10.0 Myr. 

The major conclusions of our initial paper were that 
the instability efficiently limits Mars’ mass while essen- 
tially setting its geologic accretion timescale without 
disturbing Earth and Venus’ formation, in reasonable 
agreement with geochemical studies arguing that the 
planet finished forming within the first few Myr after 
nebular dispersal (Dauphas & Pourmand 2011; Tang & 
Dauphas 2014; Kruijer et al. 2017). Our most success- 
ful outcomes occurred in simulations that most closely 
matched Jupiter’s modern eccentricity, and those where 
the instability was delayed 1-10 Myr. With shorter de- 
lays, the truncated terrestrial disk tended to spread back 
out and produce under-mass Earth and Venus analogs. 
However, none of our simulation sets provided good 
matches to the inner planets’ low degree of orbital exci- 
tation (an outstanding problem in most formation mod- 
els: Raymond et al. 2009; Lykawka & Ito 2019) and Mer- 
cury’s low mass (although generating such planets likely 
requires modifications to the C&W98 disk conditions: 
Chambers 2001; Lykawka & Ito 2017, see further dis- 
cussion in § 2.5). 

In a follow-on study (Clement et al. 2019b, hereafter 
Paper 2), we essentially reran the simulations from Pa- 
per 1 utilizing a code that accounts for imperfect ac- 
cretion by generating collisional fragments (Chambers 
2013, we use this same code in our current study, see 8 
3.1). We leveraged the same N&M12 instability models, 
however we deviated from the methodology of our ini- 
tial paper by stopping and discarding simulations where 
Jupiter and Saturn exceeded Ps/P; = 2.8. In addi- 
tion to the classic C&W98 disk, we also tested a nar- 
row annulus of embryos and planetesimals confined be- 
tween 0.7 and 1.0 au (Hansen 2009, hereafter referred to 
as the H09 disk) in both instability and control (giant 
planets on static orbits) models. The major conclusions 
of Paper 2 were that collisional fragmentation aids in 
lowering the eccentricities and inclinations of Earth and 
Venus analogs, and that the H09 annulus is compatible 
with the early instability scenario. However, while a few 
systems displayed levels of dynamical excitation com- 
parable to that of the real terrestrial system, such out- 
comes were rare. Similarly, while small Mercury analogs 
formed as the result of fragmenting collisions in some of 


1 While the original Nice Model (Tsiganis et al. 2005; Gomes 
et al. 2005; Morbidelli et al. 2005) only considered the four known 
outer planets, recent modifications include an additional one or 
two primordial ice giants to increase the probability of a simula- 
tion finishing with the correct number of planets, and reduce the 
amount of time powerful secular resonances spend in the asteroid 
belt and inner solar system (Brasser et al. 2009) by forcing Jupiter 
and Saturn's semi-major axes to evolve in step-wise manner upon 
the ejection of the additional planet(s) (Nesvorny 2011). 


our models (e.g.: Asphaug & Reufer 2014), the orbits 
of these planets were systematically too close to Venus’ 
(see Clement et al. 2019a). 


2.3. Terrestrial disks derived from high-resolution 
embryo formation models 


Initial conditions similar to our C&W98 and H09 
disk are typically justified in the literature by the re- 
sults of dust evolution models (Birnstiel et al. 2012) 
and semi-analytic studies of runaway planetesimal ac- 
cretion (Wetherill & Stewart 1989; Kokubo & Ida 1996, 
1998; Chambers 2006). Our investigations in Paper 1 
and Paper 2 indicated that the instability’s efficiency of 
limiting Mars’ mass without disturbing the other plan- 
ets’ formation and orbits is related to the partitioning 
of mass between embryos and planetesimals in each re- 
gion (Mgm»/Mpin). To better ascertain the authen- 
tic values of Mgmp/Mpin at various regions of the ter- 
restrial disk around the time of nebular dissipation, we 
performed high-resolution planetesimal accretion simu- 
lations in Clement et al. (2020a, hereafter referred to 
as the C20 disk). Our N-body models began with r ~ 
100 km planeteimsals, utilized a GPU-accelerated inte- 
gration package (Grimm & Stadel 2014) and included 
algorithms designed to mimic the effects of the decaying 
gas disk (Morishima et al. 2010). The main findings of 
our study were that embryos in the Earth/Venus region 
grow to ~0.3-0.4 Mg, and most of the small planetes- 
imals in the vicinity are accreted within the life of the 
gas disk. Contrarily, embryos do not form in the asteroid 
belt given the slow accretion timescales. These results 
are largely consistent with other recent high-resolution 
modeling efforts (Carter et al. 2015; Walsh & Levison 
2019; Woo et al. 2021a), and we validated the compat- 
ibility of the inferred disk structures within the early 
instability scenario in Clement et al. (2021b, hereafter 
Paper 3). In that work we found that the dominant ef- 
fect of the updated initial conditions was to shorten the 
planets! geologic growth timescales. In this manuscript, 
we utilize initial conditions based off these C20 disks in 
all of our new simulations. 


2.4. Updated instability evolutions 


Our analyses in Paper 1 and Paper 2 demonstrated 
how the reduction of Mars’ mass is related to the proper 
excitation of Jupiter's eccentricity; in particular the 
magnitude of its’ fifth eccentric mode es; (the term 
related to the solar system’s gs eigenfrequency: No- 
bili et al. 1989; Morbidelli et al. 2009). However, the 
adequate replication of this quality is a low-likelihood 
event in simulations of the N&M12 instability that do 
not over-excite Saturn's eccentricity (Nesvorny & Mor- 
bidelli 2012; Deienno et al. 2017). While classic stud- 
ies of the Nice Model relied on the findings of one- 
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dimensional, fixed-viscosity hydrodynamical disk mod- 
els that Jupiter and Saturn most likely to be captured 
in a 3:2 MMR (Masset & Snellgrove 2001; Morbidelli 
et al. 2007; Pierens & Nelson 2008), recent work demon- 
strated a broader spectrum of potential evolutionary 
pathways (including capture in the 2:1 with primordial 
eccentricity excitation: Pierens et al. 2014). In Clement 
et al. (2021d, hereafter referred to as the C21d insta- 
bility model) we performed a broad investigation of the 
untested parameter space applicable to the eccentric, 
primordial 2:1 Jupiter Saturn resonance and found many 
cases that greatly increase the probability of replicating 
many properties of the outer solar system (especially 
Jupiter and Saturn's precise eccentricities) when com- 
pared to the N&M12 model. Thus, we concluded that 
the scenario represents a viable evolutionary path for 
the solar system. Around half of the terrestrial planet 
formation simulations we performed in Paper 3 used the 
C21d instability model. While N&M12 also investigated 
the 2:1 Jupiter-Saturn resonance, they did not consider 
the possibility that the gas giants’ primordial eccentric- 
ities were elevated.. While the differences between the 
3:2 and 2:1 cases were mostly minor in terms of their 
effects on the inner solar system, we did note a boosted 
efficiency of Mercury-analog formation in our 2:1 cases. 
Most of these analogs were embryos liberated from the 
terrestrial disk during the instability. 


2.5. Dynamical avenues for Mercury’s origin 


Similar to many notable statistical studies of terres- 
trial planet formation leveraging N-body simulations 
(e.g.: Raymond et al. 2009; Jacobson & Morbidelli 2014; 
Izidoro et al. 2015), we largely neglected Mercury’s for- 
mation in Paper 1, Paper 2 and Paper 3. As discussed 
above, forming Mercury analogs with the correct mass 
and orbital offset from Venus requires the use of spe- 
cific initial conditions in the inner disk (a S 0.7 au, 
e.g.: Chambers 2001; Lykawka & Ito 2017, 2019), exotic 
planetary migration schemes during the gas disk phase 
(e.g.: Batygin & Laughlin 2015; Raymond et al. 2016; 
Broz et al. 2021) or low-probability collisional geome- 
tries (e.g.: Benz et al. 2007; Asphaug & Reufer 2014; 
Jackson et al. 2018; Chau et al. 2018). Studying Mer- 
cury’s formation within the early instability scenario is 
the primary goal of this paper. 

As a starting point for our current investigation, we 
use the results of a pair of recent companion investiga- 
tions by our group designed to identify viable forma- 
tion avenues for Mercury. In Clement et al. (2021a, 
hereafter C21a) we found that well-spaced (~30-40 mu- 
tual Hill Radii: Ry) systems of 4-6 ~Mars-mass short- 
period planets are easily destabilized in a manner that 
often leaves behind a relic Mercury analog with the ap- 
propriate mass, orbit and composition. In Clement & 
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Chambers (2021, hereafter C&C21) we identified a nar- 
row subset of mass-depleted inner disk parameters (in 
terms of total mass, Mgmp/Mpin and surface density 
profile) that boosts the in-situ production efficiency of 
Mercury-like planets. While we did not initially provide 
a physical motivation for these structures, in subsequent 
work (Clement et al. 2021c) we argued that both could 
be the product of proto-Earth and Venus first forming 
closer to the Sun before migrating out near the end of 
the gas disk's life. 

'The disparate interpretations of various geochemical 
and observational constraints in the literature could be 
invoked to either support or conflict with our proposed 
Earth/Venus outward migration scenario (for a com- 
plete discussion of these caveats we refer the reader to 
Clement et al. 2021c). As the precise structure of the 
primordial solar nebula remains rather unconstrained, it 
is our intention to avoid presenting a thorough justifica- 
tion of the likelihood of our chosen initial conditions oc- 
currence within the larger context of the young solar sys- 
tem's global evolution. Instead, our strategy here is to 
begin our study with the inner disk configurations that 
are seemingly most likely to produce adequate Mercury- 
analog planets. In this manner, the question of whether 
these proposed scenarios can be disentangled from the 
processes of terrestrial planet formation and the giant 
planet instability remain largely unaddressed. For in- 
stance, the C21a models mostly considered fully-formed 
versions of Earth, Venus and Mars, and only presented 
a small number of simulations that incorporated a H09 
annulus of embryos and planetesimals. The subsequent 
sections present the results of new simulations of the 
early instability scenario that modify the C20 disks uti- 
lized in Paper 3 with inner disk components derived from 
C21a and C&C21, and test both the N& M12 and C21d 
instability models. 


3. METHODS 
3.1. New simulations 


We utilize the same general methodology for generat- 
ing initial conditions, triggering the instability and se- 
lecting the evolutions that best replicate the modern 
outer solar system as described in Paper 3. Simply put, 
our full "instability plus terrestrial planet formation" 
simulations are launched on a number of compute cores, 
and simulations are continuously restarted with a new 
set of initial conditions if our algorithm detects that the 
post-instability state of the outer solar system is unsat- 
isfactory as determined by Ps/P; and ej. Specifically, 
we discard simulations that do not excite e; to greater 
than 0.03, and those that exceed Ps/P; = 2.5. Simi- 
larly, we stop any simulation where an ej output is less 
than 0.01 (half of Jupiter's modern minimum eccentric- 


ity) at any point after the instability. Integrations that 
obtain a successful instability are run up to t = 200 
Myr and saved for analysis. We also perform an addi- 
tional 10 Myr simulation of the final planetary system 
(without small bodies) using a higher coordinate out- 
put frequency to calculate the secular frequencies and 
magnitudes via frequency modulated Fourier transform 
(Sidlichovsky & Nesvorny 1996). 

Discarding all simulations that exceed a Jupiter- 
Saturn period ratio of 2.5 severely limits our sample 
size of evolved systems when compared to our previous 
studies (Paper 1; Paper 2; Paper 3) that included sim- 
ulations that finished in the 2.5 < Ps/Pj < 2.8 range. 
Indeed, of the thousands of systems we initialized for 
this work, only around <10% adequately excite Jupiter's 
eccentricity without pushing Saturn past its 5:2 MMR 
with Jupiter in the immediate aftermath of the insta- 
bility. Nearly three quarters of that subset of successful 
runs subsequently exceed Ps /P 7 = 2.5 via residual mi- 
gration or follow-on dynamical instabilities. However, 
through this process we are still able to obtain a sam- 
ple of approximately 20 fully evolved systems for each 
set of initial conditions tested. This allows us to make 
comparisons between the distributions of outcomes in 
the different simulation batches with a modest degree of 
statistical significance. For this reason though, it is im- 
portant to acknowledge that the strongest conclusions 
from our modeling work are those derived from compar- 
isons of the different disk (i.e: grouping all H09 annu- 
lus runs and scrutinizing them against the results of all 
C&W98 disks) or those contrasting the disparate insta- 
bility models. 

It is valid to criticize our methodology of combining 
the results from a suite of instability realizations (none 
of which exactly replicate the outer solar systems' true 
early evolution) as being idealized. However, alterna- 
tive methodologies that utilize interpolation algorithms 
to force the giant planets to follow ideal evolutions that 
are hand-selected from large suites of instability simula- 
tions (Nesvorny et al. 2013, 2014a,b; Roig & Nesvorny 
2015; Roig et al. 2016; Deienno et al. 2016; Nesvorny 
et al. 2021; DeSouza et al. 2021) are plagued by the same 
inherent weakness. It is impossible to know whether 
the selected instability simulation represents the true 
evolution of the giant planets or not. In our approach 
(here, as well as in Paper 1; Paper 2; Paper 3), we at- 
tempt to select a range of instabilities that most closely 
replicate the modern dynamical configuration of Jupiter 
and Saturn. While Uranus and Neptune's dynamical 
perturbations in the inner solar system during its for- 
mation (Chambers 2001; Levison & Agnor 2003), and 
today (Nobili et al. 1989; Laskar 1990), are minor; ex- 
citation of planetesimals in the terrestrial disk is closely 
related to Jupiter's free eccentricity (Raymond et al. 


2009), and the Jupiter-Saturn period ratio consequen- 
tially sets the location of secular resonances that scatter 
and remove material from the asteroid belt (Minton & 
Malhotra 2010; Walsh & Morbidelli 2011; Clement et al. 
2020b). By combining a number of different evolutions 
that replicate these qualities in broad strokes, our work 
aims to present a large sample of potential early evolu- 
tions of the young solar system; each of which possess 
a disparate sequence of ice giant encounters, peak gas 
giant eccentricities, ice giant low perihelia passages, Jo- 
vian semi-major axis jumps, and residual migration dis- 
tances. However, it will be important to study whether 
any of our most successful instabilities (in terms of in- 
ner solar system constraints) present systematic issues 
for small body populations in the outer solar system 
(e.g.: Nesvorny et al. 2013, 2014a; Deienno et al. 2014; 
Nesvorny 2015) in future work. 

Table 1 summarizes the important assumptions of our 
different simulation sets, along with the parameters uti- 
lized in a number of simulation sets from our past early 
instability papers used here for the purpose of compari- 
son. All simulations leverage a version of the Mercury6 
Hybrid integrator (Chambers 1999) that includes modi- 
fications (Chambers 2013) designed to model imperfect 
accretion and collisional fragmentation (Genda et al. 
2012; Leinhardt & Stewart 2012; Stewart & Leinhardt 
2012). When the angle and velocity of an impact place it 
into the erosive regime, the fragmenting mass is divided 
into a number of equal-mass particles with m > 0.005 
Me that are ejected at uniformly-spaced directions in 
the collisional plane with v œ 1.05ve5-. While the inclu- 
sion of collisional fragmentation does not appreciably 
affect the broad orbital architecture and mass distribu- 
tion of the resulting planetary systems (e.g.: Deienno 
et al. 2019), we include it in our simulations primarily for 
the purposes of tracking the evolution of Mercury's core 
mass fraction (CMF). The integration timestep in our 
simulations is 3 days, objects are considered to merge 
with the Sun at 0.05 au, ejections occur at Q — 100.0 
au, and we include additional forces to account for gen- 
eral relativity (Saha & Tremaine 1992). All terrestrial 
objects begin with densities of 3.0 g/cm?. Each individ- 
ual computation considers four separate regions of the 
solar system (see figure 1) derived from previous simu- 
lations: 


e (1) Inner (a < 0.7 au) terrestrial disk. 
'This section of the terrestrial region is alterna- 
tively referred to as the *Mercury-forming region" 
throughout our manuscript. All of our simula- 
tions utilize one of two sets of initial conditions 
determined to be successful in our previous inves- 
tigations of Mercury's formation (8 2.5). Simula- 
tions using the C21a model distribute five, 0.05 
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Ma embryos between 0.3-0.7 au with circular, co- 
planar orbits and randomly assigned semi-major 
axes that separate them from each other and ob- 
jects in the outer terrestrial forming disk by ~30- 
40 Ry. In successful realizations, the instability 
destabilizes this system in a manner such that a 
single Mercury analog remains as the lone sur- 
vivor. Contrarily, runs using the C&C21 initial 
conditions distribute 20 equal-mass embryos and 
200 equal-mass planetesimals between 0.3-0.7 au. 
The total mass of planet-forming material is set to 
0.5 Mg in all models, and equally distributed be- 
tween the embryo and planetesimal populations. 
Thus, inner disk embryos are ten times more mas- 
sive than the planetesimals. Eccentricities and in- 
clinations are randomly selected from near-circular 
Rayleigh distributions, and semi-major axes are 
assigned in a manner such that the surface density 
profile of the region increases linearly as X œ r. 


(2) Outer (0.7 « a « 4.0 au) terrestrial disk. 
The initial conditions for the Venus, Earth and 
Mars-forming regions used in our new simulations 
are identical to those employed in Paper 3 with 
the notable exception that, here, we remove ob- 
jects with a « 0.7 au and replace them with the 
new inner disks described above. The orbits and 
masses of our outer disk embryos (m > 0.01 Ma) 
are the same in all of our simulations, and are 
taken directly from the GPU-accelerated C20 sim- 
ulations of runaway growth at t — 8 Myr (5 Myr 
after nebular dispersal). Conversely, each time a 
simulation restarts we randomly select a new dis- 
tribution of 861 planetesimals by drawing from the 
8 Myr C20 population of planetesimal (m « 0.01 
Ma) orbits with a > 0.7 au. The objects are as- 
signed equal-masses such that the cumulative mass 
of all planetesimals in different radial bins of the 
outer disk is the same as that of the original GPU 
models. 


(3) Giant planet resonant chain. As in Paper 
1, Paper 2 and Paper 3, we pre-evolve our giant 
planet resonant chains in the presence of a plan- 
etesimal disk (the primordial Kuiper Belt) up to 
the point of a close-encounter (d < 3.0 Ry) be- 
tween two giant planets before inputing our ter- 
restrial disks (t;,,,, in table 1) and evolving the 
entire solar system through the instability period. 
Through this process, we generate a large sam- 
ple of outer solar systems on the verge of insta- 
bility (typically the gas giants semi-major axes 
diverge within 10-100 kyr after restarting these 
simulations). The integration timestep for these 
pre-instability simulations is 50 days. As in Pa- 


per 3, we test both a five-planet, 3:2,3:2,3:2 chain 
(N&M12 model) and a five-planet, 2:1,4:3,3:2,3:2 
C21d model where the gas giants’ initial eccentric- 
ities are ~0.05. Each time a compute core restarts 
an instability with a new planetesimal distribu- 
tion, it also randomly draws a new set of outer 
solar system initial conditions from our sample. 


e (4) Primordial Kuiper Belt. The planetesimal 
disk in our initial pre-instability simulations con- 
sists of 1,000, equal-mass objects with anep+1.5 < 
a « 30 au. The surface density of the disk falls off 
as r+, and the total disk mass is set to 35.0 Mg 
in our N& M12 instabilities, and 20.0 Mẹ in our 
C21d model. In most cases, the disk is depleted 
in mass by around a factor of two by the time we 
input our pre-evolved outer solar system into our 
terrestrial planet formation simulation. Thus, our 
full simulations consider a solar system composed 
of ~1,500-1,750 objects. 


3.2. Constraints 


Our current analyses combine a number of success 
criteria developed and used in our previous studies 
of terrestrial planet formation. We also employ two 
constraints introduced in Nesvorny et al. (2021) that 
eliminate certain redundancies and degeneracies related 
to commonly-used metrics such as angular momentum 
deficit and radial mass concentration (e.g.: Laskar 1997; 
Chambers 2001). Table 2 summarizes the 11 impor- 
tant criteria that are relevant to our current study of 
the early instability. However, we remind the reader 
that our outer solar solar system constraints (e55 and 
Ps/ Pj) are necessarily satisfied by all simulations as a 
result of our computational pipeline (8 3.1). 

To analyze a given system, we must first determine 
which simulation-generated planet (or planets) is an 
analog of each of the actual terrestrial planets. We be- 
gin by selecting the two most massive terrestrial planets 
in each system, and denoting them as Earth and Venus 
analogs. In contrast to our previous work, we do not an- 
alyze simulations that only form a single planet (these 
infrequent cases are discarded and not included in the 
discussion of our results). 

Our first three criteria scrutinize Mercury and its dy- 
namical relationship to Venus, and are similar to the 
metrics employed in C&C21. As simulations occasion- 
ally form multiple small planets in the Mercury region, 
rather than picking one and scrutinizing its mass ratio 
with Venus, we are more interested in the total mass of 
all planets interior to Venus (Mue/My = Y; Milai < 
ay )/My). We consider a system to be successful if it fin- 
ishes with Mme/My < 0.2. Systems where Mercury’s 
mass is greater than 3.6 times that of the real planet 


(0.2 Ma, i.e.: those with a single Mercury and a mas- 
sive Venus), as well as those where no Mercury analog 
originated as an embryo are not included as success- 
ful. While our upper limit on Mercury's mass is rather 
large, we find it necessary to provide adequate statis- 
tics given the small number of Mercury analogs finish- 
ing with masses less than twice that of the actual planet 
in our study. As a typical Venus analog has a mass of 
~0.5-1.0 Ma in our models, and each system is initial- 
ized with a cumulative planetesimal and embryo mass of 
0.25-0.5 Ma in the nominal Mercury-forming region, to 
be successful a system must still lose a significant frac- 
tion of the material in the region (typically to merger 
with the Sun or Venus). 

Another distinctive feature of Mercury and its rela- 
tionship to Venus is the innermost planets! dynamically 
isolated orbit (explained in detail in C21a). Indeed, even 
the lowest mass Mercury analogs reported in the ter- 
restrial planet formation literature tend to possess or- 
bits that are far too close to Venus’ (e.g.: Chambers 
2001; Lykawka & Ito 2017, 2019). To asses our simu- 
lations' capacity for replicating this peculiar solar sys- 
tem quality, we require the orbital period ratio between 
Venus and the largest Mercury analog be greater than 
2.25. While Mercury's moderately excited inclination 
is also noteworthy for being larger than that of any of 
the other planets in the solar system, this quality is not 
particularly challenging to match in numerical simula- 
tions (see discussion in: Roig et al. 2016; Clement et al. 
2019a; Nesvorny et al. 2021). Finally, Mercury's dis- 
tinctive massive iron core (70-8096 of its total mass as 
estimated by the MESSENGER mission: Hauck ot al. 
2013) presents an interesting constraint for our Mercury 
analogs. While it is possible that physical or chemical 
processes altered the compositions of Mercury's precur- 
sor planetesimals (e.g.: Ebel & Alexander 2011; Wurm 
et al. 2013; Kruss & Wurm 2018; Johansen & Dorn 
2022), we are keenly interested in understanding the de- 
gree to which fragmenting impacts in our simulations 
might alter Mercury's CMF during the giant impact 
phase (Benz et al. 1988; Asphaug & Reufer 2014; Ebel & 
Stewart 2017). To be considered successful, a Mercury 
analog (the most massive planet interior to Venus) must 
finish with a CMF greater than 0.5. Our methodology 
for tracking core and mantle transfers in fragmenting 
collisions is described in detail in Clement et al. (20192). 
Simply put, we assume that collisional fragments are 
first produced from the mantle of the projectile. If the 
impacting body's mantle is totally eroded, fragments 
are then produced from the projectile's core, followed 
by the mantle of the target, and finally its core. For all 
simulation sets, we assume that all embryos and plan- 
etesimals have CMF — 0.3 around the time of nebular 
dispersal (we refer to this as “time zero" throughout the 


Set Gin-Gout (au) Nemb Npin Memb tot (Me) Moin,tot (Me) tinsto (Myr) Neim Frag 
Comparison models from past work 
C&W98/CJS (Paper 1) 0.5-4.0 100 1000 2.5 2.5 N/A 50 
C&W98/CJS (Paper 2) 0.5-4.0 100 1000 2.5 2.5 N/A 100 N 
H09/CJS (Paper 2) 0.7-1.0 400 0 2.0 0.0 N/A 100 Y 
Comparison instability simulations from past work 
C&W98/N&M12 (Paper 1) 0.5-4.0 100 1000 2.5 2.5 1 18 N 
C&W98/N&M12 (Paper 2) 0.5-4.0 100 1000 2.5 2.5 1 22 Y 
H09/N&M12 (Paper 2) 0.7-1.0 400 0 2.0 0.0 10 17 Y 
C20/C21d (Paper 3) 0.48-4.0 23 954 2.25 2.33 5 16 N 
New instability simulations in this paper 
C20+C&C21/N&M12 0.3-4.0 26 1061 1.69 2.22 5 16 Y 
C20--C21a/N&M12 0.3-4.0 30 861 1.82 2.09 5 20 Y 
C20--C&C21/C21d 0.3-4.0 28 1061 1.69 2.22 5 15 
C20+C21a/C21d 0.3-4.0 48 861 1.82 2.09 5 22 Y 


Table 1. Summary of initial conditions for complete sets of terrestrial planet formation simulations. The columns are as follows: 
(1) the terrestrial (outer+inner) disk combination and giant planet instability model used for each simulation set (recall that 
N&M12 refers to the five planet instability model with Jupiter and Saturn in the 3:2 resonance, and C21d refers to the five planet 
model with the gas giants in a 2:1 MMR, see the text and figure 1 for definitions of the different disks), (2) the inner and outer 
edges of the terrestrial forming disk, (3-4) the total number of embryos and planetesimals, (5-6) the total mass of the embryo 
and planetesimal components, (7) the instability timing in Myr (ie.: the time giant planets are added to the simulation), (8) 
the total number of integrations comprising the set, and (9) whether or not (Y/N) the computations were performed using the 
collisional-fragmentation version of Mercury (Chambers 2013). The acronyms CJS and EJS refer to Jupiter and Saturn being 
placed on Circular and Eccentric (modern) orbits, respectively. 


Criterion Actual Value Accepted Value Source Notes 
Mercury Constraints 
Mme/Mv 0.0674 <0.2; >0.01 C&C21 crit. B Calculated as $5 Mi(a; < av)/Mv 
Py /Pme 2.55 22.25; qv > Qu C&C21 crit. C Pme of most massive planet with a < av 
CM Fme 0.7-0.8 > 0.5 C&C21 crit. D See Paper 2 for CMF calculation 
Other Terrestrial Constraints 
(e; i)gv 0.027 <0.055 Nesvorny et al. (2021) (e,t)ev = (ev + eg + siniv 4- sinig)/4 
Aagv 0.28 au < 0.4 Nesvorny et al. (2021) AaEv = QE — av 
TO 30-100 Myr >30 Myr Paper 1 crit. C See Kleine et al. (2009) 
TMa 1-10 Myr «10 Myr Paper 1 crit. B See Dauphas & Pourmand (2011) 
Mya/Mg 0.107 «0.3; » 0.01 Paper 1 crit. A Calculated as X` Mi(a; > ag)/Mg 
Man,f/MAn,o < 0.001 « 0.02 Clement et al. (2019c) 
Outer Solar System Constraints 
€55 0.044 > 0.022 Q tinse +10 Myr N&M12 crit. C 
Ps/ Pj 2.49 2.3-2.5 N&MI2 crit. D « 100 kyr between 2.1-2.3 


Table 2. Summary of the success criteria used in our analyses. The rows are: (1) the ratio of the total planetary mass (m > 
0.01 Ma) remaining inside of Venus’ orbit to Venus’ mass, (2) the Mercury-Venus orbital period ratio for system's with adequate 
values of Mme/Mv, (3) Mercury's final core mass fraction in systems forming analogs, (4) the total eccentricity and inclination 
excitation of the system's Earth and Venus analogs, (5) the Earth-Venus semi-major axis spacing, (6-7) the time for Earth and 
Mars to accrete 9096 of their mass, (8) the ratio of the total planetary mass remaining outside of Earth's orbit to Earths' mass, 
(9) the fraction of asteroidal mass remaining after the instability simulation, (10) the magnitude of Jupiter's fifth eccentric mode 
and (11) the Jupiter-Saturn period ratio. 
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remainder of the text). 

We use the new constraints introduced in Nesvorny 
et al. (2021) as a starting point for constraining our 
Venus-Earth-Mars systems. The authors of that paper 
rightly recognized that the standard metrics employed 
throughout the literature, namely angular momentum 
deficit (AMD: Laskar 1997) and radial mass concentra- 
tion (RMC: Chambers 2001), possess certain degenera- 
cies that can potentially complicate a rigorous statisti- 
cal interpretation of simulation outputs. We analyze 
these degeneracies in greater detail, and provide fur- 
ther justification for the updated Nesvorny et al. (2021) 
metrics in a forthcoming paper (Clement et al. 2023). 
To briefly summarize these arguments, while AMD and 
RMC scrutinize the orbital excitation and spacing of 
the inner planets as a group, the actual solar system is 
more accurately described by the contrasting properties 
of the larger and smaller terrestrial worlds. In partic- 
ular, Earth and Venus’ orbits are remarkably dynami- 
cally cold, while those of Mercury and Mars are not. We 
scrutinize this property directly by averaging the eccen- 
tricities and inclinations of the two most massive planets 
(Earth and Venus analogs) as: 


ev +ep +siniy +sinig (1) 
4 

To be successful, a system must attain a value of (e, i) gy 

within a factor of two of the solar system statistic. Sim- 

ilarly, we measure the Earth-Venus orbital spacing as: 


(e, i) BV = 


Aagy = ap — ay (2) 


Our remaining constraints on Earth and Mars analogs 
(Te, TMars and Mya/Mg; table 2) are similar to those 
employed in our past studies of the early instability sce- 
nario. Here, we define the most massive planet beyond 
Earth's orbit as the system's Mars analog when cal- 
culating TMars, and combine all planets in the region 
when determining Mma/Mpg. If no planets form beyond 
Earth, the system fails both analyses. More specifically, 
we require Earth analogs attain >90% of their total mass 
in >30 Myr (consistent with isotopic dating of the tim- 
ing of the Moon-forming impact: Yin et al. 2002; Wood 
& Halliday 2005; Touboul et al. 2007; Kleine et al. 2009; 
Rudge et al. 2010; Zube et al. 2019), and the largest 
planet exterior to Earth (Mars analog) accrete 9096 of 
its mass within the first 10 Myr of the simulation (based 
off constraints on the planets' accretion derived from 
analyses of the Martian meteorites: Dauphas & Pour- 
mand 2011; Tang & Dauphas 2014; Kruijer et al. 2017; 
Costa et al. 2020). However, it is important to note 
that the precise temporal history of Earth's accretion 
and the corresponding timing of the Moon-forming gi- 
ant impact are still debated (e.g.: Barboni et al. 2017; 
Thiemens et al. 2019). As some systems form multiple 


small planets in the Mars region, in contrast to our pre- 
vious works, we require the total mass of planets in the 
Mars-region be no more than 3096 the mass of the sys- 
tem's Earth analog. As with our Mercury constraints, 
systems forming no Mars analog that originated as an 
embryo are not considered to be successful. 

It is also important to point out that the large ec- 
centricities of Mars (0.093) and Mercury (0.21) in the 
solar system are important qualities to be replicated in 
any formation model. Past work on terrestrial planet 
formation and the evolution of the solar system dur- 
ing the giant planet instability have found that these 
values are common results in simulations (Roig et al. 
2016; Clement et al. 2019a; Lykawka & Ito 2019; DeS- 
ouza et al. 2021). For this reason, we do not include spe- 
cific constraints on the eccentricities and inclinations of 
Mercury and Mars in our analyses. Indeed, the median 
eccentricity of Mercury (0.17) and Mars (0.12) analogs 
in our study is in good agreement with those of the real 
planets. 

While our simulation planetesimals are each individ- 
ually more massive than the estimated cumulative mass 
of the modern asteroid belt, we are still interested in 
quantifying the degree of depletion in the belt in our 
early instability models. Indeed, if the asteroid belt was 
heavily populated with planetesimal mass after nebular 
dispersion as assumed here (Hayashi 1981; Bitsch et al. 
2015), it must have subsequently been depleted by some 
3-4 orders of magnitude in order to match its present low 
mass (Petit et al. 2001). Though long-term dynamical 
evolution does transfer mass from the belt to the near- 
earth asteroid region, it is estimated that the total mass 
of asteroids has only been reduced by a factor of ~two 
or so over the past 4 Gyr (Minton & Malhotra 2010). 
While our simulations’ number of initial planetesimals 
in the region (575) hinders us from scrutinizing depletion 
with the same accuracy employed in other recent high- 
resolution studies (e.g.: Roig & Nesvorny 2015; Deienno 
et al. 2016, 2018; Clement et al. 2019c), and the precise 
initial mass of the belt being poorly constrained (e.g.: 
Raymond & Izidoro 2017b; Dermott et al. 2018), as a 
first-order measure of success we require our systems 
finish with fewer than 296 of the initial asteroidal mass. 
Thus, assuming an additional depletion factor of ~50% 
over the subsequent 4 Gyr, our successful simulations 
exclusively comprise instabilities capable depleting the 
belt by at least two orders of magnitude. 

An obvious way in which our classification scheme can 
break down is in the assessment of a system dominated 
by three or more massive planets where the most massive 
two (Earth and Venus analogs) do not inhabit neighbor- 
ing orbits. In this case, an overly massive Mars analog 
might be misclassified as an Earth analog, and the true 
Earth analog (i.e.: the planet or planets in between the 


two most massive bodies in the system) would not be 
considered by any of our analyses. This could poten- 
tially artificially boost our rates of success for My, /Mg 
and artificially deflate the number of systems satisfying 
our Aagy constraint. However, as the Nice Model dis- 
turbance efficiently removes material from the Mars and 
asteroid belt regions, the two most massive planets in all 
of our instability models are almost exclusively neigh- 
boring planets in the Earth and Venus-region. Indeed, 
we verified that only three of our instability systems fin- 
ished with a planet with m > 0.3 Ma in between the two 
most massive terrestrial worlds. Thus, this degeneracy 
does not appreciably affect the interpretation of the re- 
sults for our instability models. Conversely, such config- 
urations occur relatively frequently in our CJS models. 
For this reason, we restrict Earth and Venus analogs in 
CJS models to planets with a « 1.3 au. 

It is also possible that a smaller (e.g.: Mars-mass) 
planet forms between Earth and Venus, and would there- 
fore not be considered in our analyses. Upon closer in- 
spection we find that, while such outcomes do occur 
in all of our simulation sets, the additional unclassified 
planet typically inhabits an unstable orbit that would 
lead it to eventually merge with Earth or Venus. As 
we would not expect such an instability to radically 
shift Aagy or (e, i) gv (see, for example: DeSouza et al. 
2021), and we do not impose an upper limit on Earth's 
accretion timescale (75), ignoring such planets should 
not appreciably affect the interpretation of our results. 


3.3. Comparison Simulation Sets 


Throughout the subsequent sections we compare and 
contrast our results with those of a number of simulation 
sets from our past studies (see table 1). As discussed in 
§ 1, outer terrestrial disks (i.e: beyond 0.7 au) utilized 
in N-body studies of the various formation scenario can 
essentially be grouped in two separate categories: ex- 
tended disks (C& W98) and truncated (annulus) disks 
(H09). While lower order variations on these basic con- 
cepts are obviously important for a detailed study of a 
particular model (e.g.: our updated C20 extended disk 
initial conditions for our early instability-model-focused 
study), it is also important to place the specific results of 
a study such as ours in the appropriate broader context 
with clear comparisons to alternative scenarios. For this 
reason, we refer to two sets of 1109 style simulations (one 
with circular giant planets, and one utilizing a N&M12 
instability model) initially presented in Paper 2. These 
models utilized the same fragmentation code (Chambers 
2013) described in § 3.1, and were initialized with 400 
equal-mass embryos distributed between 0.7 and 1.0 au. 
As annulus models such as these are far more success- 
ful at replicating the terrestrial (e, i) gy value than ex- 
tended disk models (Nesvorny et al. 2021), understand- 
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ing whether this advantage is maintained when the Nice 
Model disturbance is accounted for will be an important 
focus of our subsequent analyses. Similarly, we also refer 
to control runs with and without fragmentation (CJS), 
and the most successful instability simulations from Pa- 
per 1, Paper 2 and Paper 3 in order to understand how 
the different disk structures and calculation methodolo- 
gies (namely our collisional fragmentation code) affect 
our results. To maximize consistency, we utilize simu- 
lations considering 1 and 5 Myr instability delays from 
those works, and only include simulations that finish 
with Ps/P; < 2.5 and ess > 0.022. In order to pro- 
vide a sufficient number of systems for statistical analy- 
sis, we combine C&\W98/N&M12 style simulations that 
were initialized with three and four primordial ice giants 
from both Paper 1 and Paper 2 (see figure 1). 


4. RESULTS 


'Table 3 summarizes the percentages of systems in each 
of our simulation sets (including both reference mod- 
els from our previous papers and new calculations per- 
formed in this work) that satisfy our various success 
metrics (8 3.2). The subsequent sections describe the 
results of our combined analyses, and are organized as 
follows: § 4.1 focuses on Mercury, § 4.2 discusses our 
final Earth-Venus systems, and § 4.3 briefly addresses 
Mars and the asteroid belt (discussed in more depth in 
our previous studies: § 2). 


4.1. Mercury Analogs 


The primary goal of this study is to assess the com- 
patibility of the new inner disk structures (figure 1) 
identified in C21a (lone survivor model) and C&C21 
(mass-depleted embryo/planetesimal component) with 
our early instability framework. Unsurprisingly, each 
simulation set including planet-forming material in the 
inner, a < 0.7 au region of the terrestrial disk performs 
remarkably better than any of our previous reference 
sets when scrutinized against our Mercury-specific con- 
straints: Mye/Myv, Py / Pye and CM Fye. Indeed, 
Mercury- Venus systems are essentially ubiquitous in the 
majority of these models. In this manner, figure 2 plots 
an example temporal evolution of a simulation that uti- 
lized our C20--C21a/N&M12 initial conditions, and si- 
multaneously satisfied a number of our success criteria. 


4.1.1. Fewer Mercury Analogs from the 2:1 Jupiter-Saturn 
resonance 

Our simulations that initialize Jupiter and Saturn in 
a 2:1 MMR with elevated eccentricities (C21d model) 
represent a notable exception to the overall trend of 
Mercury analog generation in our integrations. Specifi- 
cally, these new simulation sets possess the highest frac- 
tion of system finishing with no mass interior to Venus 
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Set Myue/My Py/PmMe CMF ue (e,i)ev Aapv To Tma Mmwa/Mg  Man,;/MAnB,o 
Comparison control models from past work 
C&W98/CJS (Paper 1) 0 12 N/A 55* 8 85 9 o** 0 
C&W98/CJS (Paper 2) 6 11 2 83* 33 80 33 2** 0 
H09/CJS (Paper 2) 36 5 2 100* 56 27 5l 55 N/A 
Comparison instability simulations from past work 
C&W98/N&M12 (Paper 1) 0 5 N/A 44 27 93 12 16 33 
C&W98/N&M12 (Paper 2) 13 0 0 50 72 90 24 40 50 
H09/N&M12 (Paper 2) 17 5 0 76 70 20 64 52 N/A 
C20/C21d (Paper 3) 12 6 N/A 37 37 35 45 18 50 
New instability simulations in this paper 
C20+C&C21/N&M12 38 100 54 46 38 46 45 38 92 
C20+C21a/N&M12 11 78 35 52 58 94 40 35 94 
C20+C&C21/C21d 26 80 40 20 13 46 22 46 100 
C20+C21a/C21d 13 75 75 36 22 77 68 22 95 


Table 3. Summary of the percentage of systems in each of our various simulation sets (table 1) that satisfy each of our success 
criteria (table 2). Values in bold highlight the most successful simulation set according to each metric. The columns are as 
follows: (1) the name of the simulation set, (2) the ratio of the total planetary mass (m > 0.01 Mẹ) remaining inside of Venus’ 
orbit to Venus’ mass, (3) the Mercury-Venus orbital period ratio for system's with adequate values of Myre/Mv, (4) Mercury’s 
final core mass fraction in systems forming analogs, (5) the total eccentricity and inclination excitation of the system's Earth and 
Venus analogs, (6) the Earth-Venus semi-major axis spacing, (7-8) the time for Earth and Mars to accrete 9096 of their mass, 
(9) the ratio of the total planetary mass remaining outside of Earth's orbit to Earths' mass, and (10) the fraction of asteroidal 
mass remaining after the instability simulation. Note that (*) CJS simulations would require a later Nice Model instability to 
excite the giant planets’ eccentricities, and this event would also significantly increase (e,i)gv (e.g.: Kaib & Chambers 2016). 
Additionally, we note that (**) the most massive planet formed in C& W98/CJS disks tends to be near Mars’ modern semi-major 
axis. Thus, our classification algorithm typically labels it as the Earth analog. To correct for this, we require Earth and Venus 
analogs in these simulations have a<1.3 au. 


(10/15 C20+C&C21 simulations and 14/22 C20--C21a 
systems). Contrarily, only 6 total simulations from both 
of our new instability batches using the N&M12 model 
entirely deplete the Mercury region. This contrasts with 
our findings in Paper 3. There, we noted an increased 
fraction of Mercury analogs formed when Jupiter and 
Saturn were trapped in the primordial 2:1 resonance. 
However, those models did not initialize any planet 
forming material interior to 0.5 au, and the few Mercury 
analogs that did form began the simulation as embryos 
in the Venus-region or, more commonly, the Mars region. 
As Jupiter's eccentricity is close to its modern value for 
the entire duration of these simulations, perturbations 
from the migration of the v5 resonance (which begins at 
low eccentricity and inclination around a ~ 0.75 au for 
P;/Ps = 2.0) are stronger throughout the duration of 
the simulation. While a detailed analysis of this result 
is beyond the scope of this paper, we remind the reader 
that reasonable Mercury analogs (according to all three 
of our metrics) are still obtained in these models (pink- 
shaded points in the subsequent figures). However, the 
majority of these planets have masses that are 2-3 times 
greater than Mercury's actual mass. Moreover, in 8 4.2 
we will argue that the dynamical structure of the Earth- 
Venus system is also potentially inconsistent with the 2:1 
instability model. Thus, the fact that these models also 


struggle to produce Mercury analogs further strengthens 

our arguments in favor of the 3:2 (N&M12) model. 

4.1.2. Matching the Mercury- Venus period ratio with C21a 
disks 

Figure 3 plots the distribution of Mercury-Venus mass 
and period ratios attained in all of our new simulations, 
compared with selected reference models from our past 
work (table 1) out to a period ratio of 5.0. To highlight 
the most successful analog systems, we plot a box that 
bounds the systems finishing with 2.25 < Py / Py. < 3.0 
and 0.034 < Mye/My < 0.2. These limits are based on 
our success criteria for the values (table 2), the Mercury- 
Venus 3:1 MMR, and half the current value of My ec /My. 
Given the number of outer solar system constraints al- 
ready applied to these simulations prior to analysis, it 
is encouraging that there is at least one reasonable Mer- 
cury analog in this box for each instability model (3:2 
and 2:1) and inner disk structure (C21a, C&C21 and 
H09 annulus). 

While some of our Mercury-Venus systems exceed 
Py /Pme = 5.0 (discussed below), we exclude these from 
this figure as they obviously represent poor outcomes. 
'Though it is clear from the distribution of plotted val- 
ues that the outcomes of our new instability simulations 
(blue and pink shaded points) span the complete range 
of possible Py /Pme values, they also tend to finish with 
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Figure 2. Semi-Major Axis/Eccentricity plot depicting the evolution of a successful system in the C20--C21a batch. The size 
of each point corresponds to the mass of the particle (because Jupiter and Saturn are hundreds of times more massive than the 
terrestrial planets, we use separate mass scales for the inner and outer planets). For reference, the embryo initial embryos that 
go on to become Venus, Earth and Mars are color-coded gold, blue and red, respectively. Collisional fragments are annotated 
as pink triangles. The final terrestrial planet masses are Mme = 0.13 (Mme,ss = 0.055) My = 1.10 (My,ss = 0.815), Mg = 
0.85 and Mma = 0.15 (Mma,ss = 0.107) Ma. The additional important success metrics satisfied by this simulation include: 
Py /Pme = 2.49, Aagy = 0.40 and Tẹ = 50.5 Myr and Mya /Mg = 0.30. However, Earth and Venus’ orbits are overly excited 
((e,i)Egv = 0.21: a common problem in our study, see section 4.2.2), Mars’ formation timescale is slightly too long (Tma = 17 
Myr) and the asteroid belt is not sufficiently depleted (MAB,f/MAB,o = 0.040). 


Mure/My ~ 0.1-0.3; slightly greater than the solar sys- 
tem value of 0.067. While this is clearly a marked im- 
provement from all of the extended disk (C& W98) mod- 
els from Paper 1 and Paper 2, it is interesting that the 
H09 (annulus) simulations more consistently yield Mer- 
cury analogs with the appropriate mass. However, with 
the exception of one instability model (light green point 
on top of the red star in figure 3; our best Mercury- 
Venus system), no annulus model finishes with Py / Pure 
in excess of the solar system value. 


To further expound on the challenges involved in repli- 
cating the Mercury-Venus period ratio in our models, 
we plot the cumulative distribution of all final system 
Py /Pye statistics for all of our instability simulation 
batches in figure 4. 'To bolster statistics in this plot, we 
include simulations that finish with 2.5 < P;/Ps < 2.8 
that are not considered in any of our other analyses. It 
is again clear from this figure that our reference models 
from Paper 1 and Paper 2 that truncate the inner ter- 
restrial disk around the vicinity of Venus' modern orbit 
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Figure 3. Final values of Mme/Mv versus the Mercury-Venus period ratio, Py / Py. for all Mercury-analogs formed in our 
present study (see explanation in $ 3.2) compared with those from our past studies. Simulations that do not include a giant 
planet instability (CJS) are annotated with squares. Instabilities where Jupiter and Saturn begin in a 3:2 MMR (N&M12 model) 
are plotted in shades of blue, while shades of red or pink triangles indicate 2:1 (C21d) models. All annulus (H09) simulations 
are plotted in shades of green. The red star indicates the solar system value. The black square highlights the region bounded 
between a period ratio of 2.25 (our success criterion) and 3.0 (not to exceed the 3:1 MMR), and a mass ratio of 0.034 (half the 
current value to ensure stranded planetesimals are not included) and 0.2 (our success criterion). 


struggle to produce Mercury analogs in general. More- 
over, the rare analogs themselves are systematically too 
close to Venus. Furthermore, a second trend emerges 
upon closer inspection of the differences between our 
two inner disk models (C21a and C&C21). While the 
majority of our C&C21 simulations that place a col- 
lection of 20 embryos and 200 planetesimals in the in- 
ner disk finish with Mercury- Venus period ratios signif- 
icantly in excess of the actual value, the distribution of 
results for our C21a systems cluster more tightly around 
Py /Pme = 2.6 (solar system value) in figure 4. Despite 
our C& C21 simulation sets possessing slightly improved 
success rates for the first three metrics listed in table 3, 
for this reason, we assess our C20+C21a initial condi- 
tions to be the most successful in terms of their ability 
to generate the authentic Mercury-Venus system. 

In our initial study of the “lone survivor" scenario 
where ~4-6 Mars-mass embryos interior to Venus desta- 
bilize and leave behind a single Mercury analog though a 
sequence of fragmenting collisions (C21a) we exclusively 


modeled the giant planets on their current orbits. A 
major shortcoming of those simulations was a tendency 
of the final Mercury-Venus analogs to have excessively 
large orbital period ratios, however we speculated that 
this might be resolved if the event transpired during gi- 
ant planet migration. Our new simulations confirm this 
suspicion. Specifically, 13/20 total simulations in our 
C20--C21a/N&M12 set (lone survivor model with 3:2 
version of the Nice Model) form Mercury analogs with 
2.0 < Py/Pme < 3.0. While only one of these analogs 
simultaneously satisfies our relatively strict constraint 
on Mercury's mass, it is clear from figure 3 that very 
few systems are successful in this manner generally (we 
note that there is still a large amount of uncertainty 
in collisional fragmentation models that might also ac- 
count for some of Mercury's mass depletion, see § 4.1.3 
for additional discussion). Moreover, an additional two 
of these 13 systems finish with two Mercury analogs, nei- 
ther of which is more massive than 0.2My. Thus, a late 
collision with Venus would be the only event separating 
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Figure 4. Cumulative distribution of final Py / Pare values for Mercury-Venus systems formed in giant planet instability models. 
The black dashed vertical line indicates the solar system value. Instabilities where Jupiter and Saturn begin in a 3:2 MMR 
(N&M12 model) are plotted in shades of blue, while shades of red or pink indicate 2:1 (C21d) models. All annulus (H09) 


simulations are plotted in shades of green. 


these simulations from success. 

'To better understand how our new instability simula- 
tions provide better matches to Py /Pme than our orig- 
inal C21a models, we inspect the frequency at which 
each specific proto-planet survives to become the sys- 
tem's Mercury analog, as well as the average time of, 
and reason for the loss of the other four embryos. In 
our former simulations that only considered the evo- 
lution of 5 proto-planets interior to the seven other 
planets on their modern orbits, the innermost embryo 
survived as the system's final Mercury analog 44% of 
the time. Contrarily, this was only the case in 32% of 
our C20--C21a/N&MI12 instability simulations. In fact, 
we also observe instances where a rouge planetesimal 
or embryo from the Venus-forming region is scattered 
inward onto a Mercury-like orbit. Additionally, with- 
out perturbations from the instability included, it is far 
more likely for all proto-planets to combine into a single, 
overly massive Mercury analog. In C21a we observed no 
cases where the innermost embryo merged with Venus 
or Earth, and only a single instance in over 200 simula- 
tions where the second proto-planet merged with one of 
the fully-formed terrestrial planets. Contrarily, 9 of the 
embryos closest to the Sun in our 20 instability simula- 
tions combine with Venus; thus reducing the probability 
of forming a Mercury analog with an excessively large 


value of Py/Pye. Thus, we argue that the incorpora- 
tion of an instability model improves the likelihood that 
Mercury forms at the correct semi-major axis by remov- 
ing excessive embryos with low semi-major axes. 


4.1.3. Enhanced C M Fme in instability simulations 


Figure 5 plots the cumulative distribution of Mercury 
CMFs for all of our fragmentation simulations (refer- 
ence models from Paper 2 and new runs from this work). 
Notably, regardless of the specific structure of the ter- 
restrial disk, instability simulations produce a greater 
fraction of high-CMF Mercury alongs than the corre- 
sponding CJS simulations. This is unsurprising given 
that the instability-induced eccentricity excitation of the 
terrestrial disk provides high-speed collisions that have 
the potential to fall in the fragmentation regime. Cou- 
pled with the aforementioned improved Py / Py. results 
for certain simulation sets, this result illustrates how an 
early instability might provide a compelling explanation 
for Mercury's peculiar orbit and internal structure. 

It is also clear from figure 5 that the inferred value 
of CM Fme (70.7) lies outside the range of values pro- 
duced in any of our H09 annulus models (both with and 
without giant planet migration). This is largely a con- 
sequence of the fact that Mercury essentially forms as 
a stranded embryo in these models. Once Mercury is 
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Figure 5. Cumulative distribution of final CMFs of Mer- 
cury analogs formed in simulation sets utilizing collisional 
fragmentation (Chambers 2013). The black dashed vertical 
line indicates Mercury's actual approximate inferred CMF 
(Hauck et al. 2013). Instabilities where Jupiter and Saturn 
begin in a 3:2 MMR (N&M12 model) are plotted in shades of 
blue, while shades of red or pink indicate 2:1 (C21d) models. 
All annulus (H09) simulations are plotted in shades of green. 


ejected from the annulus region, its accretion is essen- 
tially over. Thus, these Mercury analogs tend to have 
less overall time and correspondingly fewer accretion 
events to alter their CMFs. While it is certainly possible 
that the Mercury's precursor planetesimals and embryos 
were altered in CMF via chemical or physical processes 
prior to its ejection from the region, it would be diffi- 
cult to explain how this process did not also affect the 
Earth or its precursors since both planets must neces- 
sarily originate in the same annulus. However, as the 
nature and size of Venus’ core remains largely uncon- 
strained (e.g.: Aitta 2012; O'Neill 2021), future explo- 
ration (e.g: DAVINCI+ and VERITAS) will undoubt- 
edly be key in providing improved constraints for our 
terrestrial formation models. 

In general, similar fractions of Mercury analogs in our 
new instability simulations (rightmost blue and pink 
lines in figure 5) attain CMFs in excess of 0.5. It is 
important to note that the percentages provided in the 
fourth column of table 3 report the fraction of all simula- 
tions in the respective set that finish with a high-CMF 
Mercury analog. Thus, sets of initial conditions that 
struggle to form such planets in general (specifically our 
2:1 C21d instabilities) finish with lower scores, even if 
the distribution of CMFs in figure 5 is similar to those 
of the more successful batches. 

While the CMFs of Mercury analogs in our C21a simu- 
lations that only consider five proto-planets in the inner 
disk are typically altered in a series of fragmenting inter- 
actions between a pair of specific embryos, the CMF evo- 


lution of inner disk objects in our C&C21 disks is often 
exceedingly complex. Figure 6 plots one example evolu- 
tion for a system from our C20+C&C21/N&M12 batch 
that satisfies all but one (Tma) of our success criteria. In 
total, 39 collisional fragments (color-coded pink in the 
top panel) are ejected from this analog over the dura- 
tion of its growth. One intriguing aspect of this analog's 
evolution is the fact that, prior to being permanently 
enhanced in CMF, the proto-Mercury is involved as the 
smaller object in five high-velocity, erosive glancing col- 
lisions and re-ejected as a mantle-only fragment. As the 
original projectile in such interactions is considered the 
"first fragment" by our algorithm, it is always assigned 
the initial mantle portion of the ejected material. Simi- 
larly, several near equal-mass accretion events involving 
other fragments composed of mostly core-material radi- 
cally alter the analog's CMF in the positive direction. 

It is important to consider how the random assign- 
ment of core and mantle material to the particles pro- 
duced in fragmenting collisions when post-processing 
our results can artificially enhance or reduce the final 
CMF of a specific Mercury analog in a specific simula- 
tion. If a hypothetical final surviving, high-CMF planet 
in the inner disk originated as a core-only fragment af- 
ter a catastrophically destructive impact, it is easy to see 
how the system might have been unsuccessful if we had 
simply assigned mantle material to that specific frag- 
ment, rather than core-material. Thus, we are far more 
interested in the distribution of final CMFs in figure 5 
than we are in, for example, the degree of mutual exclu- 
sivity between high C'M Fme and other constraints in 
individual systems. 


4.2. Earth- Venus System 
4.2.1. Earth's growth and the Moon-forming impact 


As demonstrated in table 3 and our past studies of 
the early instability (8 2), prolonging Earth's accretion 
(Tẹ > 30 Myr) is not a challenging constraint to satisfy. 
While Earth grows rapidly if the planet-forming mate- 
rial is concentrated in an annulus (such models often 
argue that a five planet inner solar system was rapidly 
assembled, and then survived in a quasi-stable state un- 
til the Moon-forming impact, e.g.: Johansen et al. 2021; 
Broz et al. 2021; Izidoro et al. 2021b), accretion in our 
extended disk models is typically more prolonged. 

Figure 7 analyzes the nature of the final giant impacts 
on Earth analogs (the nominal Moon-forming event) in 
our various instability simulation sets. While previ- 
ous studies have scrutinized the connection between the 
types of impacts produced N-body studies of terrestrial 
planet formation (e.g.: Kaib & Cowan 2015; Quarles & 
Lissauer 2015; DeSouza et al. 2021; Woo et al. 2022) 
and the particular types of Moon-forming events found 
to be consistent with orbital and chemical constraints 
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Figure 6. Top panel: Perihelia and aphelia of the four terres- 
trial planets from a simulation in our C204-C&C21/N&M12 
set. The pink lines plot the q/Q of collisional fragments 
ejected from the Mercury analog (black line) that alter its 
CMF during its growth. Middle and Bottom panels: Mass 
and CMF evolution of the same Mercury analog, respectively. 


via hydrodynamical modeling (e.g.: Canup 2004b; Cuk 
& Stewart 2012; Canup 2012; Reufer et al. 2012; Lock 
& Stewart 2017; Carter et al. 2020) we did not analyze 
our simulations in this manner in our past studies that 
were more focused on Mars’ formation. 

In general, there are two types of proposed impact 
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Figure 7. The top panel plots impact parameters 


(b) and velocities (vimp/vesc) for all roughly equal-mass 
(Mrneia/Mearth > 0.7) Moon-forming impacts (defined here 
as the final collision between Earth and an embryo) in our 
new instability simulations. The solid square denotes the 
preferred parameters from the analysis of Canup (2012). The 
second panel is similar to the top panel, except here we only 
display final giant impacts where Mrneia/Meartn < 0.5. The 
solid lines mark the preferred geometry of Canup (2004b), 
and the dashed lines indicate the successful impacts from 
the rapid-rotating model of Cuk & Stewart (2012) (see text). 
'The third and fourth panels plot the cumulative distributions 
of the final giant impact times (ter) and Theia:Earth mass 
ratios in our simulation sets. Instabilities where Jupiter and 
Saturn begin in a 3:2 MMR (N&M12 model) are plotted 
in shades of blue, while shades of red or pink indicate 2:1 
(C21d) models. 
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scenarios that might have led to the formation of the 
Earth-Moon system: those where the impactor, Theia, 
has a mass approximately equal to that of Mars (the so- 
called canonical scenario, e.g.: Benz et al. 1986; Canup 
2004b; Asphaug 2014) and those preferring a collision 
between roughly equal-mass objects (e.g.: Canup 2012). 
As our C20 initial outer terrestrial disks are dominated 
by three, m ~ 0.3 Ma embryos at time zero, we naively 
expected equal-mass Moon-forming events to be com- 
mon occurrences in our simulations. However, it is clear 
from the upper two panels that this is not the case. In- 
deed, equal-mass large accretion events (defined here as 
Mrneia/Mearth > 0.7) are far more common on Venus 
early in our simulations, while Earth's accretion is typi- 
cally more prolonged, and involves the addition of mul- 
tiple Mercury-Mars mass objects (see additional discus- 
sion in Paper 3). Indeed, over ~80% of the final giant 
impacts in all of our different instability simulations have 
projectile:target mass ratios (MTheia/MEartn in the bot- 
tom panel of figure 7) less than ~0.2. 

In the canonical model for the formation of the Moon 
(e.g.: Hartmann & Davis 1975; Cameron & Ward 1976; 
Benz et al. 1986), the angular momentum of the Earth- 
Moon system is considered to be a conserved quan- 
tity. High resolution simulations of this scenario utiliz- 
ing smoothed particle hydrodynamics (SPH) codes indi- 
cate that a low-velocity (vii /vesc S; 1.05) collision at an 
oblique angle (7445?) is most consistent with constraints 
from the modern systems’ mass partitioning and total 
angular momentum (Canup & Asphaug 2001; Canup 
2004a,b). A slight variation on this model was proposed 
in Reufer et al. (2012), where the authors advocate for 
a slightly more energetic (Vimp/Vese ~ 1.2) hit-and-run 
impact that leaves behind a Moon analog predominantly 
derived from proto-Earth material. While such a sce- 
nario potentially requires Theia be more water-rich and 
possibly originate in the main belt (Jackson et al. 2018), 
given the similarities between the Reufer et al. (2012) 
and Canup (2004b) models we mostly consider them 
jointly when analyzing our simulations. It is clear from 
the distributions plotted in the second panel of figure 7 
that similar events are common in our models. However, 
it is important to interpret the rather low rate at which 
each model produces the precise impact geometry (i.e.: 
yielding a point inside the box) within the appropri- 
ate context of the highly stochastic giant impact phase. 
Moreover, while the the distribution of Moon formation 
times (taz) in the third panel of figure 7 includes many 
early events (t < 30 Myr) that are potentially inconsis- 
tent with isotopically inferred ages (e.g. Kleine et al. 
2009), impacts that are good analogs to the preferred 
canonical impact disproportionately occur later in sim- 
ulations. Indeed, the median value of tg; for 0.5 « b « 
0.85, Mrheia/ MEarth « 0.5 and Uimp/ Uesc « 1.25 is 63.1 


Myr. 

A variation of the canonical, Mars-mass impactor 
model was proposed by Cuk & Stewart (2012). The au- 
thors proposed a high-speed impact involving a rapidly 
rotating proto-Earth and an impactor that strikes at 
an orientation retrograde to Earth's spin. Furthermore, 
the model does not require the system angular momen- 
tum be conserved, and instead argues that the rapidly 
rotating proto-Earth-Moon system can be spun-down 
through the evection resonance with the Sun. The 
dashed lines in the second panel of figure 7 denote the 
preferred impact parameters from the Ćuk & Stewart 
(2012) analysis, and demonstrate feasibility of such a 
Moon-forming event occurring within the early instabil- 
ity framework. It is interesting that similarly energetic 
interactions occur with the greatest frequency in our 
C21d instability models that initialize Jupiter and Sat- 
urn in a 2:1 MMR with elevated eccentricities. This is 
most likely the consequence of an increased tendency of 
terrestrial over-excitation in these instability scenarios. 
Indeed, both C21d sets have significantly lower (e, i) gy 
success rates (20 and 36%; table 3) than their N& M12 
instability model counterparts (46 and 5296). We elab- 
orate further on these trends in the subsequent section. 

In summary, while our simulations produce reasonable 
matches to the proposed impact geometries of a num- 
ber of Moon-formation scenarios in the literature, the 
most common final giant impact on Earth analogs in 
our models most closely matches the canonical impact 
scenario of Canup (2004b). While perturbations from 
the Nice Model instability do excite orbits in the terres- 
trial region, exotic configurations such as the equal-mass 
model of Canup (2012) and the fast-spinning scenario of 
Ćuk & Stewart (2012) are not regular outcomes in our 
simulations. However, these results should be taken in 
the appropriate context given the fact that our analyses 
are limited by small number statistics. 


4.2.2. Replicating Earth and Venus’ cold orbits 


Figure 8 compares the distributions of final (e, 7) ev 
and Aagy values for Earth-Venus systems formed in all 
simulations utilizing our three primary disk configura- 
tions: C&W98, H09 and C20. It is obvious from the 
plotted distributions that the solar system values (yel- 
low stars) lie at the extreme of the range of simulation- 
generated values. While the annulus models plotted in 
the middle panel produce the greatest number of sys- 
tems with terrestrial-like low-(e,i)gy values and small 
Earth-Venus orbital spacings, the majority of these suc- 
cessful systems are produced in simulations that do not 
include a giant planet instability model. As discussed in 
the introduction, all terrestrial planet formation models 
must account for the effects of the giant planet instabil- 
ity. While the primary difference between the instabil- 
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success criteria ((e, i) gv versus Aagy) for our various sim- 
ulation sets. The top two panels compile simulations from 
our past studies of the early instability. Specifically, the top 
panel combines instability and static giant planet simulations 
performed using classic disk initial conditions (C& W98), and 
the middle panel displays similar data for annulus (H09) disk 
conditions. The bottom panel combines our most success- 
ful simulations from Paper 3 and those from our present 
study. As in figure 3, simulations that do not include a 
giant planet instability (CJS) are annotated with squares. 
Instabilities where Jupiter and Saturn begin in a 3:2 MMR 
(N&M12 model) are plotted in shades of blue, while shades 
of red or pink indicate 2:1 (C21d) models. All annulus (H09) 
simulations are plotted in shades of green. The yellow star 
indicates the solar system value, and the solid lines enclose 
systems that finish with both metrics between 0.5 and 1.5 
times the solar system values (note that these limits to not 
correspond to our success criteria utilized in tables 2 and 
3, and are instead selected to highlight the scarcity of solar 
system-like results). 
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ity and CJS models in all panels of figure 8 is hotter 
(e, i) gy distributions (Aapy is not particularly affected 
by the instability), reasonable analogs are still produced 
in simulations that include giant planet migration. 

Only one of our reference instability simulations from 
Paper 1 and Paper 2 that utilized the C& W98 disk fin- 
ishes in the black box around the yellow star in the upper 
panel of figure 8 that denotes outcomes falling between 
0.5 and 1.5 times the solar system values. While we can- 
not rule out the C& W98 disk simulations from Paper 2 
as incompatible with the early instability scenario on 
these grounds, it is clear that our instability simulations 
utilizing the H09 and our new C20 disks provide more 
compelling results. Clearly, the distribution of outcomes 
for the H09 instability simulations cluster more tightly 
around the solar system value. However, the fact that 
the C20 disks produce reasonable results as well makes 
it difficult to argue that one disk structure represents 
the authentic state of the inner solar system around the 
time of nebular gas dispersal, while the other does not. 
Indeed, our best C20 instability sets in this paper sat- 
isfy our constraints for (e,i)gy and Aagy around 50% 
of the time. as compared to ~70% for the annulus insta- 
bility runs. Although our analyses in 8 4.1 disfavor the 
H09 disks because they produce poor Py/Py. values, 
this might be improved if a concentrated annulus outer 
terrestrial disk structure were combined with one of our 
inner disk configurations (C21a or C&C21). 

Nevertheless, for our current purposes it seems rea- 
sonable to simply conclude that the H09 and C20 disks 
are both compatible with orbital constraints from the 
Earth-Venus system when they are evolved through the 
Nice Model instability. However, our work reinforces 
the fact that Earth and Venus' cold orbits are extremely 
challenging to replicate with embryo accretion models, 
and we note two important caveats on our overall find- 
ings: 


e Collisional fragmentation likely played a 
minor, albeit important role in Earth and 
Venus? dynamical evolution: It is also clear 
from figure 8 that our best results in terms of si- 
multaneously replicating both (e, i)gy and Aagy 
occur almost exclusively in simulations that in- 
clude collisional fragmentation (all points except 
the dark blue circles and black squares in the top 
panel). This is consistent with our findings in Pa- 
per 2. In that study, we concluded that added 
dynamical friction from ejected fragments tends 
to damp the eccentricities and inclination of grow- 
ing Earth and Venus analogs (see also: Chambers 
2013; Haghighipour & Maindl 2022). 


e No satisfactory results were obtained in 
C21d instability models where Jupiter and 
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Saturn inhabit a primordial 2:1 MMR. As 
discussed in § 4.1, our C21d-style instabilities tend 
to overly-excite embryos in the Venus-forming re- 
gion. This is the result of the v5 resonance being 
positioned in between Earth and Venus’ modern 
orbits around time zero, and being stronger than 
in our N&M12 instabilities as a result of the gas 
giants’ primordial eccentricity excitation. This is 
clearly evidenced by the difference between the 2:1 
(pink/red points in the bottom panel of figure 8) 
and 3:2 (blue points) (e,i)zy and Aagy values, 
and the success rates for these metrics provided in 
table 3. 


4.3. Mars and the Asteroid Belt 


Our previous studies in Paper 1, Paper 2 and Paper 3 
were heavily focused on the ability of our models to repli- 
cate Mars’ mass and rapid inferred accretion timescale 
(e.g.: Dauphas & Pourmand 2011). The success rates for 
TMa and Mya/Mg reported in table 3 for our new simu- 
lations largely confirm our overall conclusions from those 
previous papers regarding the viability of the early insta- 
bility scenario in terms of the Mars constraints. While 
the primary focus of our current work is Mercury’s for- 
mation, in this section we briefly build on the analyses 
of the early instability’s effects on Mars and the asteroid 
belt presented our previous studies. 


4.3.1. Systems forming no Mars analog 


In Paper 3 we noted that our new, GPU-grown C20 
disks had slightly higher rates of forming no Mars ana- 
log when compared to the C&\W98 disks considered in 
Paper 1 and Paper 2. This is the result of the total 
planetesimal masses in the Mars-region being smaller 
in these models. With fewer planetesimals available to 
damp the excited orbits of would-be Mars-analogs after 
the instability, the chances of losing all material from the 
region are higher. Our new models confirm this finding, 
and we note that the effect is particularly pronounced 
in our C21d instability models that place Jupiter and 
Saturn in a 2:1 resonance with moderately enhanced ec- 
centricities. Specifically, 41% of our C20+C&C21/C21d 
simulations form no Mars. By comparison, our runs uti- 
lizing the identical terrestrial disk configuration in com- 
bination with the N&M12 instability model completely 
deplete the Mars-region only 13% of the time. When 
combined with the other shortcomings identified in the 
previous sections, this does speak against the viability of 
the C21d model within the early instability framework. 
However, we refrain from ruling the model out entirely 
as a reasonable fraction of our 2:1 instability models do 
produce adequate Mars analogs (table 3). Interestingly, 
our H09 annulus instability models yield the lowest frac- 
tion of systems with no Mars analog (596) 


4.3.2. Consequences of the 2:1 Jupiter-Saturn resonance in 
the asteroid belt 

Both our N&M12 and C21d instabilities are quite suc- 
cessful at heavily depleting à primordially massive as- 
teroid belt. The success of our instabilities initializ- 
ing Jupiter and Saturn in the 2:1 resonance is poten- 
tially related to the fact that the giant planets' orbits 
are eccentric starting at time zero. In Paper 3 we did 
not analyze whether or not the 2:1 Jupiter-Saturn reso- 
nance negatively affects the asteroid belt's orbital struc- 
ture. Figure 9 compares the orbital inclination distri- 
bution of large asteroids in the actual belt (top panel), 
with the co-added distributions of surviving asteroids in 
all of our new N&M12 (middle panel) and C21d (bot- 
tom panel) simulations. With the exception of high- 
inclination asteroids in the inner belt that are expected 
to be removed during the giant planets’ residual migra- 
tion phase (Clement et al. 2020b), both modeled dis- 
tributions provide reasonable matches to the real belt. 
It is noteworthy that the N&M12 3:2 instabilities pro- 
duce a better match to the inclination distribution of 
asteroids with a « 2.5 au. Indeed, the ratio of asteroids 
above, to those below the vg secular resonance (0.48; a 
success criteria we used in Paper 1 and Paper 2) is rea- 
sonably close to that of the actual asteroid belt (0.08). 
However, the fact that our sample size of surviving as- 
teroids is relatively small makes it difficult to draw any 
significant conclusions from this plot. 

'To first order our work confirms the compatibility of 
the C21d model with the asteroid belt's orbital struc- 
ture. We also compared the eccentricity distributions, 
as well as the relative populations of objects in different 
radial bins, and found them both to reasonably resem- 
ble the belt's modern structure. However, co-adding the 
results of multiple simulations obviously introduces ad- 
ditional uncertainties to this analysis, and thus future 
work incorporating higher particle resolution (e.g.: Dei- 
enno et al. 2016, 2018; Clement et al. 2019c) will be 
required to fully validate these conclusions. 


5. DISCUSSION AND CONCLUSIONS 


In Paper 1, Paper 2 and Paper 3 we argued that an 
unusually early Nice Model (Tsiganis et al. 2005; Lev- 
ison et al. 2008; Nesvorny & Morbidelli 2012; Deienno 
et al. 2017) instability occurring within the first 1-10 
Myr after gas disk dissipation is capable of reducing 
the mass and formation timescale of Mars analogs (a 
problem in classic terrestrial accretion models: Wether- 
ill 1980b; Chambers & Wetherill 1998; Chambers 2001; 
Raymond et al. 2009). While such a timing for the giant 
planet instability is rather early when compared with 
other sequences of events considered in the literature 
(e.g.: Nesvorny et al. 2021), a myriad of recently derived 
constraints restrict the epoch of giant planet migration 
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Figure 9. Orbital inclination versus semi-major axis distri- 
bution of large objects (D = 50 km) in the modern asteroid 
belt (top panel), compared with all surviving asteroids in our 
new simulations utilizing N&M12 instability models (middle 
panel) and C21d models (bottom panel). The dashed lines 
annotate the locations of dominant mean motion and secular 
resonances in the belt. 
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to the first ~100 Myr after the solar system's birth (see 
§ 1 and Raymond & Morbidelli 2022, for a recent re- 
view). However, the physical process of nebular photo- 
evaporation might necessitate the solar system's insta- 
bility occur at an epoch as early as considered in our past 
studies (Liu et al. 2022). In this paper, we revisited the 
early instability scenario with new simulations aimed at 
modeling the planet Mercury's formation by incorporat- 
ing a shorter integration timestep and additional planet- 
forming material interior to Venus’ modern orbit. While 
the initial conditions used in our simulations are broadly 
consistent with the outcome of models investigating the 
effects of an outward migrating proto-Venus and proto- 
Earth during the gas disk phase (Clement et al. 2021c), 
in this paper we simply combine outer terrestrial disks 
that were tested against a number of constraints related 
to the outer three terrestrial planets in previous work 
(Paper 3) with inner disks that boosted the probability 
of forming adequate Mercury analogs as determined in 
a pair of recent papers (C21a; C& C21). This is an obvi- 
ous weakness of our current modeling work, and we plan 
to test more realistic embryo and planetesimal distribu- 
tions derived directly from a single gas-disk simulation 
in future work. 

'The main conclusion of our modeling efforts is that the 
early instability framework is consistent with all poten- 
tial pathways for Mercury's formation considered here. 
The majority of our analyses focus on differences in the 
distribution of simulation outcomes for systems initial- 
ized with various disk structures in both the Mercury- 
region and outer terrestrial disk, as well as two different 
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initial giant planet configurations. The primary results 
of varying each of these parameters are summarized be- 
low. 


5.1. The Mercury-forming region 


Our analysis of Mercury's formation focused on three 
different models developed in past work: (1) the ter- 
restrial planets form from a narrow annulus of material 
from which Mercury is scattered out of and forms as a 
stranded embryo (1109), (2) a collection of 4-6 ~Mars- 
mass planets form interior to Venus in a quasi-stable 
orbital configuration that is destabilized in a manner 
that yields many high-velocity, mantle-stripping colli- 
sions and leaves behind Mercury as the sole survivor 
(C212), and (3) Mercury forms directly within a mass- 
depleted distribution of embryos and planetesimals in 
the inner solar system (C&C21). While each model pro- 
duces at least one outstanding Mercury-Venus system, 
our best results occur in the C21a lone survivor mod- 
els. In successful simulations, instability-induced eccen- 
tricity excitation enhances the rate of high-velocity col- 
lisions among the embryos, and enhances the rate at 
which all proto-planets in the region merge with Venus. 
'This process reduces the probability of forming Mercury 
with to large of a Py /Pye. Contrarily, our H09 annu- 
lus models almost always yield Py / Pc values that are 
too small, and the period ratios produced in the C&C21 
models are systematically too large. Moreover, our H09 
disks do not produce the high-speed collisions necessary 
for sufficiently altering Mercury's CMF, while our early 
instability models more regularly match this constraint. 


5.2. The outer terrestrial disk 


Our collection of early instability simulations initially 
partition embryos and planetesimals in the Venus, Earth 
and Mars-regions in three different manners: (1) a con- 
centrated annulus of equal-mass embryos between 0.7- 
1.0 au (1109), (2) an extended disk of 100 embryos and 
1,000 planetesimals with Memy/Mpin = 1.0 at all radial 
locations (C&W98), and (3) consistent with the outputs 
of GPU simulations of planetesimal accretion in a de- 
caying gas disk (C20). While the majority of our anal- 
ysis focused on the generation of Earth and Venus’ or- 
bital spacing and degree of dynamical excitation, we also 
compared the potential Moon-forming impacts on Earth 
with the results of SPH modeling of the event (Canup 
2004b, 2012; Reufer et al. 2012; Cuk & Stewart 2012). 
In general, we find that both the H09 and C20 disks are 
capable of replicating the authentic Earth-Venus system 
after being evolved through the Nice Model instability 
(only one marginally successful result was obtained with 
the C&W98 disk). While the H09 models produce the 
best statistical results, the solar system result lies near 
the low excitation and low orbital spacing extreme of 
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the distribution of results for all modeled cases that in- 
clude an instability. Additionally, it is important to note 
that Earth’s accretion is often too rapid in the annu- 
lus models (e.g.: Kleine et al. 2009). While we observe 
Moon-forming impact geometries that are analogous to 
those preferred in a range of studies in the literature, the 
most common events most closely resemble the canoni- 
cal, low-velocity Mars-mass impactor model (Benz et al. 
1986; Canup & Asphaug 2001; Canup 2004a,b). 


5.3. Instability model: 8:2 versus 2:1 Jupiter-Saturn 
resonance 


Similar to our methodology in Paper 3, our new sim- 
ulations considered two giant planet instability mod- 
els: (1) a model where Jupiter and Saturn emerge from 
the gas disk locked in a 3:2 MMR on circular orbits 
(N&M12), and (2) a second set of computations where 
the gas giants are captured in a 2:1 resonance with el- 
evated eccentricities (C21d). While the differences be- 
tween the two models are negligible in most of our analy- 
ses (e.g.: the asteroid belt’s orbital structure, Mercury’s 
CMF, and Mars! mass and formation time), we prefer 
the N&M12 model for two main reasons. First, the 
strong resonant sweeping in the Mercury-forming region 
that is characteristic of our C21d simulations tends to 
entirely deplete the region of material. Thus, the rate of 
forming no Mercury analog is much higher with the 2:1 
Juptier-Saturn resonance. This is also the case for the 
rate of forming no Mars analog. Second, these enhanced 
perturbations also have the effect of over-exciting the 
Earth-Venus system. As a result, the C21d Earth-Venus 
orbital spacings and levels of dynamical excitation are 


systematically larger than those of both the N&M12 in- 
stability models, and the actual solar system. 

In conclusion, an early (t ~ 1-10 Myr after gas dis- 
sipation) Nice Model instability provides a compelling 
explanation for the dynamical configuration of all plan- 
ets in the solar system. While such a scenario is not the 
only plausible explanation for many of these qualities 
(e.g.: Walsh et al. 2011; Bromley & Kenyon 2017; Broz 
et al. 2021; Johansen et al. 2021; Izidoro et al. 2021b), 
all models must account for the giant planets! acquisi- 
tion of their modern orbits. Moreover, as demonstrated 
here, incorporating Mercury's formation in these types 
of models can yield additional insights and constraints. 
Future studies in the mold of our current work should 
strive to increasingly incorporate novel cosmochemical, 
geophysical and dynamical constraints to break degen- 
eracies between different formation scenarios and further 
pin down the authentic structure of the terrestrial planet 
forming disk in the solar system. 
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